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Abstract 

In this paper the isogeometric Nystrom method is presented. It’s outstanding features are: it allows 
the analysis of domains described by many different geometrical mapping methods in computer 
aided geometric design and it requires only pointwise function evaluations just like isogeomet¬ 
ric collocation methods. The analysis of the computational domain is carried out by means of 
boundary integral equations, therefor only the boundary representation is required. The method is 
thoroughly integrated into the isogeometric framework. For example, the regularization of the aris¬ 
ing singular integrals performed with local correction as well as the interpolation of the pointwise 
existing results are carried out by means of Bezier elements. 

The presented isogeometric Nystrom method is applied to practical problems solved by the 
Laplace and the Lame-Navier equation. Numerical tests show higher order convergence in two 
and three dimensions. It is concluded that the presented approach provides a simple and flexi¬ 
ble alternative to currently used methods for solving boundary integral equations, but has some 
limitations. 

Keywords: Nystrom Method, Isogeometric Analysis, Collocation, Local Refinement, Boundary 
Integral Equation, 


1. Introduction 

Isogeometric analysis has gained an enormous attention in the last decade. It offers a seam¬ 
less link between the geometrical model to the numerical simulation with no need for meshing 
in. The key idea is to use the functions describing the geometry in computer aided geometric 
design (CAGD) also for the analysis. Since CAGD models are based on boundary representa¬ 
tions, a natural integration with simulation is possible using boundary integral equations (BIE). 
Most implementations of BIEs discretize the integral equations using basis functions and numeri¬ 
cal quadrature to evaluate the bi-linear forms. In this context, the isogeometric boundary element 
method (BEM) gained much attention recently. In three dimensions, the method has been applied 
to the Laplace flU, Stokes 0, Lame-Navier SSI, Maxwell |(6j and Helmholtz [7] equations with 
different discretization methods, i.e. Galerkin [0 and collocation methods lf5l . 
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The Nystrom method, originally introduced in is an alternative for the numerical solution 
of BIEs. A unique feature of the method is that the boundary integrals are evaluated directly by 
means of numerical quadrature without formulating an approximation of the unknown fields. In 
fact, it is based on pointwise evaluations of the fundamental solutions on the boundary of the com¬ 
putational domain. The pointwise nature makes it similar to isogeometric collocation IfTOll that 
appears to be a very efficient method to solve partial differential equations [UTlI . The Ny strom 
method has been applied to potential liT2l and electro-magnetic problems lfl3l[T4| . the Helmholtz 
equation [15, J6, [ IT], Stokes flow llT8ll . to the analysis of edge cracks llT9l and elastic wave scat¬ 
tering ll20l and generally, to parabolic BIEs lUTTl . 

The implementation requires the treatment of singular fundamental solutions for which sev¬ 
eral techniques have been developed. Many of them are based on singularity subtraction [1221 
or product integration Ifl2l . Error analysis for these methods is available [f23l and the interested 
reader is referred to |]24j| or textbooks on integral equations such as llT2l or [f25l . In the presented 
implementation, the locally corrected Nystrom method lfl5ll is used for regularization. It is based 
on the local construction of special quadrature rules for singular functions li26l using composite 
quadrature rules. An overview of the method can be found in E71l and Il28ll . Mathematical proofs 
for the solvability and convergence of the locally corrected Nystrom method have been provided 
in J29) and lf30ll . 

Higher order convergence rates can be achieved with the Nystrom method based on the order of 
the chosen quadrature. A requirement for higher order convergence is the continuity of the bound¬ 
ary representation which has to be smooth [fT2l . As a consequence, standard triangulations are 
insufficient for curved geometries due to the potential kink between patches along their common 
edge. With the help of re-parametrization, patches with the required continuity may be constructed 
OTTl but the procedure is complex. This is why the Nystrom method has so far been mostly ap¬ 
plied to analytical surfaces. Its application to CAGD geometry descriptions has the advantage that 
continuity of the boundary can be better controlled, as already noted in [fl5ll . However, real-world 
geometries still contain comers and edges by design. In order to restore higher order convergence, 
the elements which define the integration regions are graded Ifl2l . For two-dimensional problems, 
several other approaches exist [32. 33, 34.35,361. While it is claimed there that these methods can 
be extended to three-dimensional problems as well, only few results are reported 071 . A solution 
for tackling singularities, arising at regions with mixed boundary conditions is presented in P8l 
and in Il39l . 

In the following sections, the isogeometric Nystrom method is introduced for the Laplace and 
the Lame-Navier equation. The formulation is based on CAGD boundary representations which 
are partitioned into integration regions for the composite quadrature. Therefor, the arbitrary se¬ 
lectable CAGD technology is only required to provide a valid geometrical mapping from the pa¬ 
rameter to the real space. The numerical scheme consists of point collocation on the surface of 
the computational domain. The regularization by means of local correction is performed with 
Bezier elements. To preserve higher order convergence, a priori grading of elements at corners 
and edges is realized in the parameter space. Without loss of generality, the procedure is adapted 
to boundary representations based on non-uniform rational B-splines (NURBS). For tensor prod¬ 
uct descriptions of surfaces the authors present a strategy for the local refinement of elements. For 
post-processing purposes, the pointwise existing results are interpolated over the boundary again 
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with the help of Bezier elements. 

In summary, the analysis with the proposed isogemetric Nystrom method offers the following 
main advantages to both, the IGA and the Nystrom community: 

• isogeometric collocation scheme for boundary integral equations, 

• free choice of CAGD technology taken for the boundary representation, 

• local refinement even for tensor product patches, 

• application of the Nystrom method to real-world geometries. 

The paper is organized in the following parts: In section [2] the boundary integral equation 
for heat conduction and elasticity is revisited. The principles and requirements of the locally 
corrected Nystrom method are explained. Section [3] describes the application and adaptation of 
the isogeometric framework to the method. Emphasis is given to NURBS surface descriptions 
and to the local refinement. In section [4] several numerical tests in two and three dimensions are 
presented. These results as well as the advantages and the limitations of the isogeometric Nystrom 
method are discussed in section [5} The paper closes with concluding remarks in section |6j 

2. Boundary Integral Equations 

Mathematical models described by mixed elliptic boundary value problems (BVP) are consid¬ 
ered in this section as well as their discretization with the Nystrom method. They are based on 
Laplace’s equation and on the Lame-Navier equation which describe isotropic steady heat con¬ 
duction and linear elasticity respectively. 

2.1. Boundary Value Problem 

Let Q. be the considered domain, u a generalized unknown and C an elliptic partial differential 
operator. Then the BVP is generalized to the following form: Lind u(x) so that 

Cu{x) =0 Vxefi 

Tu(x) = q(y) = g N (y) Vy eT n (1) 

Tr u(x) = q(y) = go(y) Vy G T D . 

Lor (jTj) the boundary F = dQ. is split into a Neumann and a Dirichlet part F = r,v U F d so that 
T N fl T/) = 0. The boundary trace 

Txu{x) — lim u{x) = u(y) xE£l,yET (2) 

x-j-y 

maps the primal field u(x) in the domain Q. to u(y) on the boundary F. The prescribed Neumann 
and Dirichlet data are denoted by g,\r and g D respectively. 

For the heat equation the partial differential operator is defined by the Laplacian A = V • V 


Cu{x) —kAu(x) 
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where k denotes the conductivity and u the temperature. The conormal-derivative maps tempera¬ 
ture u to heat flux q on the surface F and is defined by the normal derivative 

Tu(x) = kVu(y)-n(y) = q(y) xeCl,yeF (4) 

with n denoting the unit outward normal vector of Cl. 

For the Lame-Navier equation 

Cu(x) — (A + 2 ji) V • Vu(x) + /iV x (V x u(x)) = 0 (5) 

holds with u denoting the displacement and with the Lame-constants A and fl [f40ll . For elasticity 
the conormal derivative is 

Tu(x) = XV-u(y)n(y) + 2iiVu(y)-n(y) + iin(y) x (Vx«(j/)) = q(y) xeCl,yeT (6) 
which maps displacements u to surface traction q. 

2.2. Boundary Integral Equations 

The variational form of the BVP (jT|) can be solved by means of boundary integral equations 
mmm. For the solution of many physical problems, Fredholm integral equations of the first 
kind 


u(x) — (V(j>)(x) \/x e r (7) 

or the second kind 

u(x) = ((C + JC)\if)(x) Vx E r (8) 

are used. In these equations 

(Vq)(x) = j'[)(x,y)q(y)dsy Vx,y £ T (9) 

r 

denotes the single layer and 

(/Cm) (a;) = j Tll(x,y)u(y)ds y \/x,y £Y (10) 

r 

the double layer boundary integral operator. The kernel function U is the fundamental solution for 
the underlying problem which depends on the Euclidean distance r — \x — y\. In case of Laplace’s 
equation it is 

(H) 
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for QeM 3 . In case of the Lame-Navier equation it is Kelvin’s fundamental solution [[4 ll 


U ij(x,y) 


1 ^ +4 ( , X + 3y \ 

Sttr + 2/i) A+/i l ') 


( 12 ) 


represented as a tensor with i,j= (1,..., 3}. The fundamental solutions tend to infinity if r —y 0. 
Therefore the integrals in <[9j) and ( |T0| ) are and need to be regularized or evaluated analytically ll42l . 
For ( |T0| ) the regularization results in an integral free jump term 

(Cu)(x) with C — (13) 


on smooth surfaces. 

For indirect formulations with first (jTJ) or second kind integral equations ([8]), the unknowns 
(j) and y are usually non-physical quantities and only intermediate results in order to evaluate 
quantities in the interior x e £1 by means of the representation formulas 


u(x) = (V(j>)(x) and u(x) — (JC\j/)(x) \/xe£l. (14) 

Working with physical quantities u and q on the boundary F only, the BVP ([I]) is solved by means 
of a direct boundary integral formulation 

((C + JC)u)(x) — (Vq)(x) Vxer (15) 


instead. 


2.3. Discretization with the Nystrom Method 

In his original paper, Nystrom []9l| proposed the discretization of second kind integral equations 
([8]) by means of a numerical quadrature 

n 

((^l + IC)Y)(xi) ^cY(xi)+Y,Tli{xi : yj)Y{yj)wj x t e r, yj e r. ( 16 ) 

7=1 


In this equation, yj are the evaluation points of the n-point numerical quadrature and wj are their 
corresponding weights. In order to set up a linear system of equations, the quadrature sum ( p~6] > is 
collocated at distinct points x t with i — { \.... ,n \ resulting in the matrix equation 


(il + K)v< = u 

on smooth surfaces. For Fredholm integral equations of the first kind, the system is 


(17) 


V0 =u. 


(18) 
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For a direct formulation with mixed boundary conditions, a block system of equations 

X G To : f\ dD ~^Dn\ Aj/A _ f^DD —VDn\ /^g£>\ /in', 

xETy: \Vnd —K^n) \ u n J \J^nd —Vnn) \8n) ■ 

is taken like the formulation presented in ff43ll . If integral free jump terms are present, they are 
already integrated in the system matrices K. If the surface T is smooth and the kernel function U 
is regular, entries of the system matrix only consist of pointwise evaluations 

V[/,;] = U(xi,yj)wj and K[i,j\ = SijC + TU{xi,yj)wj. (20) 

For the considered applications the kernel functions are singular with c — \ on smooth surfaces. 
Moreover, the fundamental solution is undefined if x t = y } so that special treatment is necessary 
to evaluate the corresponding system matrix entries. 

Based on a technique for the construction of quadrature rules with arbitrary order for given 
singular functions presented in lf26l . the authors of llT5l developed the locally corrected Nystrom 
method for the solution of the Helmholtz equation. This particular regularization technique is 
taken for the framework presented in this paper. The main idea is to replace the contribution of 
the original kernel function in the neighborhood Q. x of the collocation point Xj with a corrected 
regular one, so that the new kernel function is defined by 


U *(xi,yj) = 


L {xi,yj) if Xie& x 

U (xi,yj) otherwise. 


( 21 ) 


The locally corrected kernel L for the collocation point Xi is computed at n corresponding field 
points yj G by solving the linear system 


£ N i(yj) L ( x i,yj)wj= J U(xj,y)Ni(y)dsy 

;=1 rntlx 


with 


i = {!,..., m}. 


( 22 ) 


Equation ( |22| ) introduces a space of m testfunctions Nj. hence the singularity of the original kernel 
is treated in a weak sense on the right hand side. The choice of the test functions in the presented 
application is discussed in section [3} Because of the singularity of U, treating the moments on the 


right hand side of fl22| ) requires regularization. This equation constructs a numerical quadrature, 
where the weights wj are not explicitly calculated but collected together with the corrected kernel 
to Wj = L(x[. yj)wj. Finally, the linear system in matrix form is 


Nw = g 


with 


N G w G M'\ g G R' n . 


(23) 


The matrix N consists of evaluations of the test-functions Ni and g contains accurately evaluated 


singular moments. Equation ( ]23j ) is numerically solved for w with LU -decomposition if n — m. 
The number m of chosen test functions iVj may be smaller or larger than the number n of sample 
points yj. In that case, a valid solution is found by means of least-squares or a minimum norm 
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solution If26ll . 

In this paper, the order of a quadrature p is defined as the degree of the highest polynomial that 
it does integrate exactly. Therefor, the numerical solution of equations ( fT7] ), ( fi~8| ) or converge 
with p + 1 = n + 1 for an n-point quadrature. Although composite Gauss quadrature is used in the 
presented formulation, the order of convergence for the method does not reach 2n + 1 due to the 
local correction Il26l . 

In practical applications, the considered surface F of the domain is approximated by F h con¬ 
sisting of non-overlapping patches r so that 


L 

T h = 1J T(24) 

l =l 

As a requirement for the locally corrected Nystrom method, a composite quadrature based on F/, is 
chosen for the discretization of the integral equation. Additionally, the quadrature is chosen to be 
of open type, which means that no quadrature points are located at the boundary of the integration 
region. This is because such points can be located at physical edges, where the integral kernel 
may be undefined or diverging. This particular choice comes with an advantage: In contrast to 
the BEM and due to its pointwise nature, the Nystrom method does not require or impose any 
connectivities between patches. Hence, non-conforming meshes are supported inherently. 

To converge with respect to the order of the underlying quadrature, the Nystrom method for 
solving integral equations of the kind ([7]), ([8]) or ([15]) requires a smooth surface. Usually, this is not 
feasible by means of a standard triangulation of F. Therefor, the authors propose the application 
of CAGD surface descriptions, which fulfill this requirement. Moreover, F/, = F which means that 
the discretization introduces no geometry error. 

3. Isogeometric Framework 

In this section the isogeometric paradigm is combined with the locally corrected Nystrom 
method. The term Cauchy data is taken to refer to quantities on the boundary F appearing in the 
discrete boundary integral equations ( |T7] ), ( fT8| ) and ( [T9] ). 

The key aspect of the isogeometric concept is to utilize the boundary representation of design 
models directly in the analysis. Thus, it has been applied to a variety of models with different 
surface descriptions such as subdivision surfaces [|44|l , tensor product surfaces JU and T-spline 
surfaces a. Since the most commonly used CAGD technology in engineering design are tensor 
product surfaces and based on non-rational B-splines (NURBS), the paper focuses on geometry 
descriptions based on this approach. However, the implementation of the Nystrom method to other 
surface representations is straight forward, since the approximation of the Cauchy data as well as 
the partitioning of the elements for the integration is independent of the geometrical parametriza- 
tion. In fact, the only requirement is a valid geometrical mapping X{r) from local coordinates 
r = (ri,. .., rj_ i) T to global coordinates x = (jci, ... ,.Vrf) T in the d-dimensional Cartesian system 
W 1 . To be precise the Gram’s determinant has to be non-singular. The corresponding Gram matrix 
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is given by 


G(r) := i\{r)i x (r) ^(d-\)x{d-\) 


(25) 


where the Jacobi-matrix 

Jx( r ) ■= with i = {l,---,d},j = {l,...,d-l} (26) 

results form the geometrical mapping X(r) ll42l . 

3.1. Geometry Representation 

CAGD objects are usually defined by a set of boundary curves or surfaces. They are generally 
referred to as patches for the remainder of this paper. A significant property is that the continuity 
within patches can be controlled by the associated basis functions. Hence, such patches are able 
to represent a smooth geometry without any artificial comers or edges. 


3.1.1. Basis functions 

The knot vector E is the fundamental element for the construction of the basis functions. It is 
characterized as a non-decreasing sequence of coordinates r,- < r- l+ \ which defines the parametric 
space of a patch. The coordinates itself are called knots and the half-open interval [r,-. r l+ \) is 
called knot span. Knot values may not be unique which is referred to as the multiplicity of a 
knot, which is then larger than one. Together with a corresponding polynomial degree p a set 
of piecewise polynomial basis functions B ip , so called B-splines are defined recursively. After 
introducing piecewise constant functions (p = 0) 




1 if r f < r < r i+i 
0 otherwise, 


(27) 


higher degree B-splines are constructed as a strictly convex combination of basis functions of the 
previous degree 


r—r: 


B,. n (r) = —— ^_i(r)+ ri+p+l ' £ m , p _i(r). 


— r 


o,p\ 


' i+p 


'i+p+l r i +1 


(28) 


Further, their first derivatives are also a linear combination of B-splines of the previous degree 


dB i , P ( r ) 


= B'Jr) = 


hP 


r i+p 


VtW 


r i+p+l — r i+l 


Bi+l, P -l{r). 


(29) 


The support supp{B/ iP } = {r,...., r i+p+ \} is local and entirely defined by p + 2 knots. B i p is 
described by a polynomial segment within each non-zero knot span (r s < r s+ 1 ) of its support. 
The continuity between adjacent segments is C p ~ m where m denotes the multiplicity of the joint 
knot. Consequently, the continuity of the basis functions B\ p at their knots is determined by the 
corresponding knot vector. Figure [Tj illustrates this relation for two quadratic B-splines. Note 
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that the continuity between the polynomial segments decreases at the double knot. In general, the 



Figure 1: Examples of two quadratic B-splines (p = 2) with different continuities at their middle knots. The knot 
vector of the left B-spline is defined by E = {1,2,3,4}, E = {1,2.5,2.5,4} defines the knot vector of the right one. 
The resulting polynomial segments are indicated by dashed lines, solid lines represents the B-splines. 

knots are arbitrarily distributed. This is emphasized by referring to the knot vector as non-uniform. 
The term open knot vector indicates that the first and last knot are C°-continuous. A special knot 
sequence is 


S = {r 0 = ••• = r p ,r p+ i = ••• = r 2p+ 1 } (30) 

where the multiplicity of all knots is equal to the polynomial order, i.e. p + 1. The resulting basis 
functions Bj p are classical pth-degree Bernstein polynomials which extend over a single non-zero 
knot span. 

3.1.2. Curves 

The geometrical mapping of B-spline curves of degree p is 


X(r) 


l- 1 


:= x(r) = £B,- P (r)c* 
f—o 


(31) 


and its Jacobi-matrix is given by 


J *(r):=£fi;, p (r)c ( - (32) 

(=0 

where I is the total number of basis functions due to a knot vector £/ and the control points Cj are 
the corresponding coefficients in physical space. The resulting piecewise polynomial curve x(r) 
inherits the continuity properties of its underlying basis functions. If they span only over a single 
non-zero knot span, i.e. £/ is of form (|30j), the curve is referred to as Bezier curve. 

In order to represent models like conic sections which are based on rational functions, weights 
Wj are introduced. They are associated with the control points yielding to homogeneous coordi¬ 
nates given by 


= (w i c i ,w i )T = (c> i )iel w 


(33) 
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Applying the mapping ( |3T| ) to c 1 ’ defines a B-spline curve x h (r) in the projective space R rf+1 . 
Hence, the geometrical mapping has to be extended by a perspective mapping y with the cen¬ 
ter at the origin of W l+1 . This is depicted in Figure [ 2 ] for a circular arc in M 2 . Its homo- 


c 


Jl 



w — 1 


X 


.h 



Figure 2: Perspective mapping y of a quadratic B-spline x h {r) in homogeneous form R 3 to a circular arc x(r) in 
physical space K 2 . 

geneous representation is defined by quadratic B-splines and the control points = (2,0,2) T , 
Cj = (2,3,2) T and c\ — (0,4,4) T in the projective space R 3 . The homogeneous vector compo¬ 
nents x w — (X|...., x h d ) 1 are mapped to the physical space by 




with 


(34) 


1—0 


The resulting rational function which describes such a curve is called NURBS (non-uniform ratio¬ 
nal B-spline). The geometrical mapping with NURBS is defined by 



(35) 


where the derivatives are defined by 



and 



(36) 
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3.1.3. Tensor Product Surfaces 

B-spline surfaces are defined by tensor products of univariate basis functions which are related 
to separate knot vectors E; and Ej. The geometrical mapping is given by 

T’(r) := ®(r) = £ £ B^ Pl (ri)B J:P2 (r 2 ) c Lj (37) 

i=0j=0 


with p i and pi denoting different polynomial degrees. The corresponding Jacobian is calculated 
by substituting the basis functions by its derivatives, alternately for each parametric direction. 
Further, a surface is labeled as Bezier patch, if Ej and Ej are of form pO] ) and the derivation of 
NURBS surfaces is analogous to NURBS curves. 

In general the surface representation by tensor products is an extremely efficient technique 
compared to other geometry descriptions [|45l . 

However, the applicability of the tensor product approach is limited. Especially in the context 
of refinement which can not be performed locally. In the scope of the present work, this drawback 
does not limit the ability for local refinement with the Nystrom method. 


3.2. Approximation of Cauchy Data 
3.2.1. Element-wise Discretization 

In order to evaluate the boundary integral properly, the patch z is subdivided into a set of 
elements f, so that 


/ 

t = (J fj. (38) 

i= 1 

A distinguishing feature of the Nystrom method is that the element-wise discretization of the 
Cauchy data is directly expressed through points defined by the applied quadrature rule on the 
geometry, rather than control points. The quadrature points are distributed with respect to their 
coordinates E, e [— 1,1] in the reference element. Consequently, the p-refinement strategy is deter¬ 
mined by the increase of the quadrature order and thus, the quadrature points on the element. To 
resolve their location y on z in physical space, the mappings A7 (f) for the reference element and 
X (r) for the geometry are consecutively applied. For a curved element z = z these mappings are 
depicted in Figure [3] 

In contrast to other isogeometric methods, the Cauchy data is not expressed in terms of vari¬ 
ables at control points c. For the analysis with the isogeometric Nystrom method, Cauchy data are 
taken from evaluations in the quadrature points y directly, without any transformation. 


3.2.2. Patch-wise Discretization 

Generally, the arrangement of f over T is independent of the geometry. But its representation 
has to be smooth within each f. As indicated in section 3.1 the continuity of the geometry is 
directly linked to the knot vector E. It is convenient, to define z by means of an artificial knot 
vector A so that E C A. To be precise, the purpose of A is not to construct basis functions but 
to organize the global partition of z properly. Consequently, a denote knots of A. Based on the 
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X&) 



Figure 3: The distribution of the quadrature points y in the reference space the parameter space r of the geometry 
and the physical space x £ R 2 . The quadrature points are indicated by circles. Squares denote control points c. 


initial discretization, i.e. E = A, the approximation quality of the Cauchy data can be improved by 
inserting additional unique knots dj into a knot span s of A, so that a s < a, < a s+ \. This procedure 
defines an h-refinement strategy that is performed in parametric space and preserves the continuity 
requirements with respect to the geometry. As indicated in Figure [4] this is independent of the 
geometry representation. 

In order to retain higher order convergence on domains with mixed boundary conditions or 
corners, the grading of elements towards such geometric locations may become necessary Ifl2l . 
Corners are easily identified in the knot vector A multiplicity of knots m = p. The grading is 
performed by subdivision of the adjacent knot spans s — \a s .a s+ \) into m elements. This is simply 
performed in the parameter space by means of knot insertion in A. The resulting n knot values cq 
with i = {1,..., n} are defined by 

q (iY 

and a f = a s+ i - (a s+ 1 - a s ) I - J (39) 

for grading towards a s or a s+ \ respectively. The exponent q is defined by 

q > P qJ> ~ with 0 < y < 1 (40) 

7 

where p q is the order of the quadrature rule and / denotes the Holder constant [H2Tl . 

In the current implementation, mixed boundary conditions are considered insofar, that either 
Dirichlet or Neumann boundary conditions along a single patch are allowed. Grading is performed 
in the vicinity of patches with different boundary conditions. 

3.2.3. Local Refinement for Tensor Product Surfaces 

The use of artificial knot vectors A is sufficient for the partition of a curve into elements f. 
However, the extension of this concept to surface representations is limited. In particular, local 
refinement is not possible if elements f are defined by a tensor product of A/ and Aj. In this 
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a/ — a s T - (ds+i &S 








X&) 





-1 1 


0 2 4 



Figure 4: The geometry of a cubic B-spline curve is defined by the knot vector £ = {0,0,0,0,2,4,4,4,4} and 
the control points c,- with i = {0,...,4}. The elements f ; with j = {1,...,4} along the curve are defined by 
A = {0,0,0,0,1,2,3,4,4,4,4} in the parameter space. Each element is equipped with a two point quadrature rule 
defined in the reference interval | £ [—1,1]. 

section, a strategy for local refinement with the isogeometric Nystrom method is explained. The 
procedure is visualized in Figure [5} 

Global refinement is performed by means of knot insertion, i.e. inserting a; in Ay as indicated 
by the dashed line in Figure [5J a). Both non-zero knot spans in A/ are subdivided. A further 
subdivision is indented to be local for each of the elements f. For the partition of an element f 
into local elements t the definition of refinement points f is adequate. Each f is defined in the 
parametric space and located either inside or on the edge of t. Inside f the refinement point defines 
the origin of a cross which is aligned to the parametric coordinate system and subdivides f into 
four local elements t. If f is located on the edge, f is subdivided into two T. Further, a local 
grid can be defined by combining several refinement points simultaneously. The described local 
refinement options are illustrated in Figure [5jb). 

In order to enable further refinement of local elements as well, all local elements are sorted in 
a hierarchical tree structure and labeled with the refinement level £. The initial refinement level 
£ — 0 refers to the global element, hence T° = f. Each node of the tree may have a different number 
/ of ancestors because of the manifold possibilities in defining f per level. The local elements if 
generated in level £ cover the complete area of the local element t^ _ 1 of the previous level 


i 



( 41 ) 
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(a) Global refinement 


(b) Local refinement £ = 1 



(c) Local refinement £ = 2 


Figure 5: Subsequent refinement of elements defined by the tensor product of A/ = {0,0,0,1,1,2,2,2} and Aj = 
{0,0,0,2,2,2}. (a) Global refinement by inserting a, = 1 into Aj. (b) Subdivision into local element by refinement 
points of the first refinement level which are denoted by circles, (c) Further local refinement by a higher refinement 
level indicated by diamonds. 


The final partition of the global element z is defined by the sum of J local elements 


= U */ 


(42) 


;=i 


related to the leafs in the hierarchical tree structure. An example of such a locally refined patch 
with two levels of refinement is depicted in Figure |5jc). The local refinement procedure involves 
the scaling and translation of the element boundaries. Details on the construction of this mapping 


due to a given set of r are found in Appendix A 


3.3. Isogeometric Nystrom Method with Local Correction 

In the presented implementation, Gauss-Legendre quadrature rules are taken. For the analysis 
in three dimensions, a tensor product quadrature is constructed as illustrated in Figure[6j However, 
it is also feasible to apply non-tensor product quadrature or numerical quadrature constructed for 
special special purposes in that context 


As mentioned in section 2.2 the integral kernels are singular if collocation and quadrature 
point coincide x t — y r Such kernels require special treatment for a correct integration. For 
the isogeometric Nystrom method, a spatial separation of quadrature points in relation to each 
collocation point is performed. An admissibility criterion 


diam (t) < rj dist (x i: yj) 


(43) 


is introduced which separates the regime with a smooth kernel function from that one with singular 
or nearly singular behavior. If ( [43] ) is fulfilled, the corresponding entries are in the scope of th e far 
field where the system matrices consist of point evaluations only. In particular, a matrix entry of 
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the discrete single layer potential ([9]) for the isogeometric Nystrom method reads 


\[i,j] = U(x i ,y j )G(rj)J x ^ j )w j . (44) 

In ( |44] ), U (xi,yj) is the evaluation of the fundamental solution with respect to the spatial coor¬ 
dinate of the collocation point x, and that of the quadrature point y t in W 1 . The evaluation of 


£2 



Figure 6: The distribution of the quadrature points y on a locally refined surface. The basis function of the geometry 
are debited by S/ = Ej = {0,0,0,2,2,2}. The partition of local elements is debited by a global insertion of a\ = 02 = 1 
and a rebnement point f = (0.5,0.5) T . 


Gram’s determinant 

G(r j) = Jda(G(rj)) (45) 

is defined by the location r j of the quadrature point yj in the parametric space. Since the quadra¬ 
ture points are given in the reference element by its coordinates — (^-j. c,/. 2 ) J , the mapping 
r j — AT (^j) to the local element t is sufficient. For the integral transformation from reference to 
local element, the Jacobian of the mapping A^ 

/(£,-) = det(jfo)) (46) 


is evaluated with respect to the reference coordinates £,• of the quadrature point y,. Finally, wj in 
(|44|) denotes the original quadrature weight. 
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The near field zone, where the integral kernels are singular or nearly singular along the affected 
local elements is defined for point combinations where ((43]) is not fulfilled. At that point, the 
matrix evaluations are performed by means of the locally corrected integral kernels as described in 


section 2.3 To perform the local correction, B-splines are taken for the polynomial test functions 
Ni = B\ p in ( [22] ) and defined on the local element t. In particular, a Bezier interpolation is proposed 
which is defined by the knot vectors of the kind 


= {-!,•• -,- 1 , 1 ,•••,!} 


(47) 


in all parametric directions. The multiplicity of the knots is chosen so that they define at least as 
many basis functions as present quadrature points on the local element. This allows the solution 


of ([ 22 ]) by means of LU -decomposition or by solving a least squares problem. 


The integrals of the right hand side in ([22]) are 


\J(xi,y)B i: p(y)ds y 


and 


Ty\J{xi,y)B Up {y)ds y . 


(48) 


If xi G t the integrals are weakly or strongly singular. In that case the single layer integral is 
subject to transformation described in PPTll and ||48l while the double layer integral is treated 
with regularization techniques presented in 


respectively. If ([43]) is not fulfilled but x, £ t, 
then the integral is nearly singular and treated with adaptive numerical integration as described 
in 0. Practically, the extent of the region where local elements are marked as nearly singular is 
determined by the admissibility factor 17 . 


3.4. Isogeometric Postprocessing 

Once the system of equations is solved, the Cauchy data exists only in the quadrature points y. 
Thus, post-processing steps are required to visualize the distribution over the whole geometry. The 
Nystrom-interpolation Ifl29 is the most accurate procedure for this task. But it requires additional 
kernel evaluations at all quadrature points, which is computationally expensive. For the isogeo¬ 
metric Nystrom method, the following approach is probably less accurate but simpler and local to 
t. Following the isogeometric concept, each element t is represented by the Bezier patch already 
constructed for the local correction. The results in each quadrature point y } are interpolated within 
each T by means of the basis functions B i p based on a knot vector of form ( |47] ). For instance, the 
primary variable in any point £ on the reference element can be calculated with 

1 

«(0 = !>,/>(£) <*■ (49) 

i =0 

In order to compute the unknown coefficients c ; the inverse of the mapping C(£)c = u is needed 
which is defined by the spline collocation matrix. Its entries are 

C \j, i] = b Up (£ 7 ) with i = {0,..., 1} and j = (0,..., J} (50) 
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where i is the number of B-spline functions and j the number of quadrature points on t. The linear 
system (|49|) can be solved directly or in a least squares sense. 


4. Numerical Results 

In this section, numerical results are provided for academic and practical problems. The results 
are critically reviewed and remarks on limitations and open topics of the isogeometric Nystrom 
method are given. 

For the numerical analysis of the convergence of the isogeometric Nystrom method, problems 
with an infinite domain Q. are solved. The fundamental solution IHx. y) with a number of source 
points x e £1 outside of Q. is applied as a boundary condition at the quadrature points y <E F. The 
problem is solved by means of Fredholm integral equations of the first kind ([7]) in order to test the 
single layer potential V and with the second kind equation ([8]) to test the double layer potential K. 
Results are given in the interior at several points x and the error is defined by 

6h = u(x) — \J(x,x) Va; e £2, x e Cl . (51) 


The relative error is 


and measured in the maximum-norm 


^rel — 


e/z 


U(a;,£c) 

e re / 1|oo. The normalized element diameter 

, f N i/(d-i) 
i _ / ^max 


(52) 


(53) 


is used for convergence plots where Aj nax is the largest length or area of all elements z and A 
the surface length or area of the whole boundary T. Hence, a step of the process denoted as 
uniform h-refinement halves the knot span of z in two dimensions (d = 2) resulting in two new 
elements. In three dimensions (d = 3), the knot spans in both parametric directions are affected 
which produces four new elements. However, the parameter h always refers to the length or area 
in M 2 or M 3 respectively. 

In the following sections, the terminus p-refinement refers to the step-wise increase of the 
quadrature order used for the simulation. For one-dimensional elements f in M 2 , a step of p- 
refinement results in the increase by one of the local quadrature points per element. For the tensor 
product quadrature being used for analysis in M 3 , this process leads to 2p—\ additional quadrature 
points on z. 


4.1. Flower 

The convergence of the method in two dimensions is tested on a smooth flower-like geometry 
solving an exterior Dirichlet problem. The whole setting is depicted in Figure [7] on the left. The 
admissibility factor for the local correction is set to rj — 2.0. In Figure [8]the results for Laplace’s 
equation with uniform / 7 -rclincmerit are shown for different orders p of the chosen quadrature rule. 
While the double layer operator K shows the expected optimal convergence of p + 1, the single 
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Figure 7: Geometries and the corresponding NURBS control grid (squares) for the test settings in two dimensions: 
flower-like boundary in infinite domain with source points x and test points x (left) and teardrop with opening angle 
a (right) 


layer V does not. In that case, the convergence reverts to a linear rate irrespective of the order of 
the integration. However, there is a significant offset depending on the quadrature rule so that a 
convergence following 

||e«/||=o = Cexp ( ^ ni) (54) 

with respect to the number of degrees of freedom n is observed for the p-refinement strategy. In the 
case of the flower-like geometry the factor s is approximately 0.55. The constant C depends on the 
size of the integration elements. This is shown in Figure [9} where h s denotes the initial, unrefined 
and dimensionless element size in the sense of equation ( |53| ). The dashed lines follow (([54])) with 
different exponential factors s For the Lame-Navier equation, the behavior is depicted in Figure [TO] 
where the Lame-parameter are defined by the elastic constants 


A = 


Ev 

(l-2v)(l + v) 


and 


E 

2(1 +V)' 


(55) 


In that example, the Young’s modulus and Poisson’s ratio are set to E — 1 GPa and v = 0.3 respec¬ 
tively. Additionally, Figure |TT] depicts the exponential convergence with respect to the degrees of 
freedom n for the ^-refinement strategy, in that case with s ~ 0.50. 

For Fredholm integral equations of the second type ([8]) solved with the Nystrom method, con¬ 
vergence theorems have been proven for the Laplace equation in two and three dimensions Ifl2l . 
The numerical results for the double layer operator in Figure [8] comply with the theory and show 
convergence rates of p + 1 or better for a quadrature rule of order p. However, for boundary inte¬ 
gral formulations in linear elasticity with the Lame-Navier equation, literature is rather sparse. In 
the straight-forward implementation described in this paper, convergence for the discrete double 
layer operator stops at a certain level of /i-rclincmcnt. 

In contrary to the Laplace equation, where the double layer operator is weakly singular, the 
boundary integral for problems in elasticity is evaluated in the sense of a Cauchy principal value 
which requires particular attention. In the presented implementation it is ensured, that the accu¬ 
racy for the local correction is beyond the measured relative error || e re /1| oo- This affects the right 
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Figure 8: Convergence rates of the indirect Laplace problem for the discrete double layer (left) and single layer 
potential (right) on the two dimensional flower like geometry. The optimal convergence rates for the lowest and 
highest order are indicated by triangles. 


hand side in equation ( f22| ). Numerical studies have shown, that a reasonable variation of the ad¬ 
missibility factor 77 does not improve the convergence significantly. An additional study, which 
has not been presented in this paper for brevity, indicated that if the admissibility factor is cho¬ 
sen to be 77 —>■ 00 and hence every integration element is corrected, the numerical results show 
optimal convergence agin. But this also means that the isogeometric Nystrom method is in fact 
replaced by the isogeometric BEM for which convergence has been demonstrated numerically in 
0. However, it is remarkable that the convergence plateau for the flower example is reached after 
exactly 3 uniform refinement steps for all quadrature orders as depicted in Figure [lOj By altering 
the element partition such that h max ~ h mui the kink has been shifted to one more refinement step, 
but still the plateau does not disappear. 

The boundary integral operator for Fredholm integral equations of the first kind is not compact 
anymore and hence, mathematical analysis is rather involved (|42|. Moreover, for problems in 
two dimensions the fundamental solution is logarithmic. The direct application of the presented 
framework to the boundary integral equation results in only linear convergence independent of the 
chosen quadrature. The performed numerical tests show this behavior for both the Faplace and 
Fame-Navier problem. For logarithmic first kind boundary integral equations on closed curves in 
M 2 it is proposed in literature to apply a transformation so that they become harmonic. Therefor, 
convergence can be restored with respect to the numerical quadrature ll23l . 

Although the convergence with / 7 -refinement for two dimensional problems is restricted, there 
is a significant offset between discretizations with increasing quadrature order. As a consequence, 
the p-refinement strategy is preferred over / 7 -rc fine merit for the analysis of practical applications 
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Figure 9: Convergence rates for /^-refinement of the indirect Laplace problem for the discrete double layer (left) and 
single layer potential (right) on the two dimensional flower like geometry. The dotted lines indicate the exponential 
function (f54|> where the aligned number is the exponential factor s 


with the isogeometric Ny strom method. Convergence follows equation (54) in that case. 
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p = i 


P = 2 


P = 3 


P = 4 


P = 5 


Figure 10: Convergence rates of the indirect Lame-Navier problem for the discrete double layer (left) and single layer 
potential (right) on the two dimensional flower like geometry. 
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+ h s —m— l/2h s —m— l/4h s 


Figure 11: Convergence rates for ^-refinement of the indirect Lame-Navier problem for the discrete double layer 
(left) and single layer potential (right) on the two dimensional flower like geometry. The dotted lines indicate the 
exponential function (|54j» where the aligned number is the exponential factor s 
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4.2. Teardrop 


For the discretization of the direct BIE ( [15] ) with the Nystrom method, singularities of the 
solution occur at domains with corners. This is, because the gradient of the solution becomes 
singular. For the isogeometric Nystrom method the problem is tackled by the grading strategy 
described in section I3T21 

To validate the algorithm, a two dimensional domain modeled like a teardrop is analyzed. 
The domain, which is shown on the right of Figure |7[ has a single corner where the opening 
angle is chosen to be a — 90°. The second kind integral equation ( fT7| ) for the exterior Laplace 
problem is solved and it’s convergence behavior investigated. In a first step, uniform ^-refinement 
is performed for different quadrature orders. Then, grading towards the comer is applied which 
is specified by equation ( [39] ) where the number of sub-elements is chosen to be m = 6 and the 
Holder constant is y = 1. To handle the small elements near the comer, the admissibility factor for 
the local correction is set to T] — 6.0. Both results are depicted in Figure 12 The convergence is 
calculated with respect to degrees of freedom n. 



P = 1 


= 2 


P = 3 


- p = 4 


-p = 5 


Figure 12: Convergence of the indirect Laplace problem on an infinite domain with one corner. The geometry is 
subject to uniform (left) and graded ^-refinement (right). Convergence rates are indicated by triangles. 

The tests on the teardrop geometry justify the grading approach. In the case of the teardrop, the 
full convergence for Fredholm integral equations of the second kind is restored. As outlined in the 
introduction, several other techniques are available to tackle such singularities. All the mentioned 
approaches may be taken into account by the presented framework for the isogeometric Nystrom 
method without any restriction. 
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4.2.1. Cantilever 

The last problem tested in two dimensions is a cantilever beam. The beam’s dimensions are 
10 m x 1 m is subject to a constant vertical top loading t y = —\ MPa. The elastic constants are E = 
29 000 MPa and v = 0.0. The solution with the isogeometric Nystrom method is compared with 
the analytic solution by Timoshenko. The vertical displacement u y is measured at the midpoint 
(10,0) of the free end of the cantilever. Grading of the integration elements towards the corners 
is realized with a Holder constant y = 1 and m = 4 sub-elements. The admissibility factor for the 
local correction is chosen to be p = 6.0. The vertical displacement for different quadrature orders 
and different degrees of freedom n are presented in Figure [13] showing excellent results for all 
orders p > 1. 



n 


Figure 13: Vertical displacement u y in [m] for the cantilever beam at its free end calculated with the isogeometric 
Nystrom method. The analysis is carried out with different quadrature orders p and degrees of freedom n. The 
analytic solution is indicated by the dashed line. 
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4.3. Torus 

A torus geometry as depicted in Figure |T4] is used to show the convergence of the boundary 
integral operators in three dimensions. The distance from the torus center to the center of the 
tube is r\ — 0.9m. The tube-radius is set to r 2 = 0.2m. Both, the Laplace and the Lame-Navier 
equation are solved. For the latter, the elastic constants are chosen to be E — 1.0 GPa and v = 0.3. 
The results are outlined in Figure [15] and Figure [T6| 



Figure 14: Toms geometry taken for the investigation of convergence in three dimensions 


For the Laplace equation, the numerical results in three dimensions comply with the theory as 
well. Similar to problems in two dimensions, optimal convergence p + 1 for the double layer oper¬ 
ator and suboptimal convergence for the single layer operator is observed. As for the convergence 
following equation (|54]), the factor is s ~ 0.30 for both, the discrete double layer operator K of the 
Laplace and the Lame-Navier equation respectively. 

However, a convergence plateau such as for the flower geometry in Figure [TO] of section 4.1 
has not been been observed. For further refinement steps the problem exceeded the number of 
degrees of freedom acceptable for such a study. The results for the three dimensional convergence 
study sketched in Figure [T6| required the application of a fast method which introduces additional 
numerical approximation. In the presented case, hierarchical matrices ("H-matrices) are taken, 
where the approximation quality of the boundary integral operators is controllable by a recursive 
criterion [50, 43] and set to at least one magnitude lower than the achieved overall accuracy. 
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Figure 15: Convergence rates for the discrete double and single layer potential on the torus for the exterior Laplace 
problem. The optimal convergence rates for the lowest and highest order are indicated by triangles. 



Figure 16: Convergence rates for the discrete double layer potential on the torus for the exterior elastic problem. 
Results are shown for the /z-refinement (left) and the ^-refinement strategy (right). On the left, the optimal convergence 
rates for the lowest and highest order are indicated by triangles. The dotted lines on the right indicate the exponential 
function (|54| where the aligned number denotes the exponential factor s 
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4.3.1. Fichera Cube 

As a benchmark for the local refinement strategy, a heat conduction problem is solved on the 
Fichera cube and compared to an FEM solution lf5lTl . The computational domain is given by 
Yl — (—1, l) 3 — [0, l] 3 and represented by 14 NURBS patches. On the Dirichlet boundary F /> the 
temperature is set to go = 0. The flux on the Neumann boundary Fy is an accumulation of known 
analytic solutions of L-shaped domains in two dimensions with respect to the xy-, yz- and the 
xz- plane. The temperature uid and the heat flux qid on the L-shape are given by 

U2d = r 2/3 sin C^-\ , with r = ^and (56) 


qid = n[i\ 


i du2d 

dti ' 


with 


i= 1,2 


(57) 


respectively. The local coordinates of the L-shape are denoted by £,• and n represents the unit 
outward normal vector. The prescribed smooth flux gw on the Fichera cube is then accumulated 
from ( |57j ) in each plane. The geometry of the problem is shown in Figure [F7] The unknown 


z 


y 




Figure 17: Fichera cube geometry, where the Dirichlet boundary is drawn red and the Neumann boundary Uv 
green. The right picture shows the corresponding L-shape in two dimensions from which the non-zero Neumann 
boundary conditions are derived. 


Flux q on Yd has a singularity at the inverted comer of the Fichera cube. In order to approximate 
q accurately towards this singular point, local refinement of integration elements as described in 
section 3.2.3 is applied. The first three refinement steps are illustrated in Figure |T8j The results of 
the isogeometric Nystrom method are compared with an /i-adaptive FEM solution, calculated by 
the Abaqus software package. The variation of the temperature u is compared along a straight line 
defined by the endpoints (0,0,1) and (—0.5, —0.5,1) on Y N . The heat flux q is compared along a 
straight line defined by (0,0,0) and (0.5,0.5,0) onFp. Both lines are drawn in Figure[l7]and the 
results are depicted in Figure |T9} The example on the Fichera cube demonstrates the ability of the 
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Figure 18: The first three step of the local /z-refinement on the Fichera cube 


isogeometric Nystrom method to tackle physical singularities by using local refinement. 
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-FEM (p = 2, /i-adaptive) -isogeometric Nystrom method (p = 3) 


Figure 19: Temperature u and heat flux q along defined paths on the Fichera Cube calculated with different analysis 
methods 
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4.3.2. Spanner 

A single-ended open-jaw spanner is considered next. The geometry is defined by a CAGD 
model which is depicted on the top of Figure[20j The elastic constants are given by £ = 2.1 x 10 5 MPa 
and v = 0.3. Homogeneous Dirichlet boundary conditions, i.e. zero displacements u — 0 are ap¬ 
plied at the open-jaw. The handle is subjected to a constant vertical loading q = (O.O.c) with 
t z — 15 MPa. This setup and the partition of integration elements resulting from the CAGD model 
is illustrated in the middle of Figure |20| 




Figure 20: CAGD model of the spanner where thin black lines delimiting the NURBS patches (top). The partition 
into integration elements and the boundary conditions are depicted in the middle, the deflection as a result of linear 
elastic analysis is sketched on the bottom. 


The spanner is analyzed with the isogeometric Nystrom method and the results post-processed 
with the strategy described in section 3.4 The deflected tool is sketched at the bottom of Figure [20} 
The vertical and horizontal deflection at the top along a straight line is shown in Figure [21] The 
results are compared with an adaptive refined FEM analysis with Abaqus. 
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— FEM (p = 2, /j-adaptive) -isogeometric Nystrom method (p = 3) 

Figure 21: Vertical deflection u z (left) and horizontal deflection u x (right) in [mm] at the top of the spanner along a 
straight line 
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5. Discussion and Remarks 


In this section, we further remark on the properties of the isogeometric Nystrom method, prac¬ 
tical considerations and on some implementation details. 

As for many other methods for solving BIEs, the accurate evaluation of singular and nearly 
singular integrals is essential. For locally corrected Nystrom methods precise singular entries are 
achieved by means of properly evaluated entries for the right hand sides of the moment equa¬ 
tion ((23]). In the presented approach, the error of the integrals in equation ( |48] ) is controllable 
and the integration is performed with an accuracy at least one magnitude lower than the lowest 
observed errors in the numerical tests of section [4] The quality of matrix entries related to nearly 
singular integration regions is controlled by the admissibility factor 77 . Depending on the com¬ 
plexity of the computational geometry, the factor is typically 17 = {2,..., 6 }. However, increasing 
7] results to a larger locally corrected region and therefor to additional numerical effort for the 
analysis. A heuristic determination criterion for 17 may assure the desired accuracy but an adap¬ 
tive integration strategy like for BEM integral kernels is not viable due to the pointwise nature of 
the method. 

The presented results for heat conduction problems and elasticity in three dimensions demon¬ 
strate the practical applicability of the method. Although tensor product surfaces are used, local 
refinement is possible. The example on the Fichera cube demonstrates the ability of the iso¬ 
geometric Nystrom method to tackle physical singularities as well. The results for real world 
examples such as the cantilever beam or the spanner are of comparable quality to the standard 
FEM. This has been outlined in section 4.3.1 and in section [4.3.2| For such practical problems, 
the p-rcfincmcnt strategy is the first choice. However, for complicated geometries this approach is 
not always suitable. Hence, the convergence behavior for Fredholm integral equations ([7]) of the 
first kind or direct BIEs like equation ( p~5| ) requires further attention to increase the robustness of 
the method. We would like to emphasize that the presented approach describes a straightforward 
implementation of the Nystrom method with a strong focus on its implementation into a isoge¬ 
ometric framework. Implementations based on kernel splitting or other techniques may perform 
better, but such formulations depend on the underlying physical problem to be solved. In three 
dimensions efficient formulations are still a matter of research [f37l ITT ]. 

Its pointwise nature greatly simplifies the implementation of the Nystrom method into fast 
summation methods. The method has been considered in the landmark paper by Rokhlin If52l for 
the fast multipole method. One of the advantages of the Nystrom method compared to other meth¬ 
ods solving boundary integral equations is, that point-wise supports are clearly related to spatial 
regions on which fast method are usually based upon. For instance, with Galerkin BEM supports 
of the basis functions may span multiple regions and therefor the partition and the treatment of 
near- and far-field system matrix entries gets more complex. The presented application makes use 
of hierarchical matrices Il53l adopting matrix operations of the HFib-library [1541 . 


6. Conclusion 

In this work, an isogeometric framework was applied to the locally corrected Nystrom method. 
The presented approach is suitable for any type of CAGD surface representation which provides 
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a valid geometric mapping from parameter to real space. Hence, the method is applicable to 
NURBS, subdivision surfaces and T-splines in a straightforward way and can be easily adapted to 
further developments and new technologies in CAGD. 

In this paper we explain the implementation of NURBS surfaces for the geometric description 
which are the commonly used in computer aided design. For such tensor product surfaces, the 
method inherently permits local /i-refinement. 

The isogeometric Nystrom method presents an attractive alternative to classical methods solv¬ 
ing boundary integral equations, its main feature being that no approximation of the unknown with 
basis functions is required. These equations are directly solved by numerical quadrature where the 
unknown parameters are located at the quadrature points. 

The isogeometric Nystrom method is characterized by a pointwise evaluation of the fundamen¬ 
tal solution on the surface. Regularization of the singular integrals is carried out by means local 
correction with B-spline basis functions. These basis functions are also used for post-processing 
purposes, where the results at quadrature points are interpolated on the surface. 

Due to its discrete, pointwise pattern, the Nystrom method is well suited for fast summation 
methods such as the fast multipole method or hierarchical matrices. 

On test examples it has been shown that the method preforms well but that its convergence 
properties are different to commonly used methods such as the boundary element method. Some 
problems that need further attention were pointed out. It is hoped that this paper gives impetus for 
further investigation in this promising alternative to classical collocation type approaches. 


32 



Appendix A. Local Element Mapping 


The mapping from the initial element to defined by a tensor product of A/ and Ay to local 
elements tg specified by refinement points fg is discussed. In a first step, each to is represented 
by means of its corner nodes, i.e. the knots of its corresponding knots span s = (,V| ..vy) 1 "- They are 
summarized in a node matrix Ao such that 

(M[s\\ A/[si + l] A/ [51 ] A f [s\ + 1]\ 

Ao = ( Aj [ 52 ] Ay[s2+1] Aj [^2 A1 ] A/[S 2 ] ■ (A.l) 

V 1 1 1 1 / 

The mapping from to to its t| j includes the translation and scaling of the corner nodes. They are 
assembled in a transformation matrix Tj,- which is defined for each t], ( - by 

( h,r\/h,r\ 0 \ 

0 h,r 2 /^0,ri tr 2 ■ (A.2) 

0 0 l) 


The last column refers to the translation t of the first corner node, whereas the diagonal entries 
are related to the lengths of the initial (Jo) and refined element (/ 1 ) in each parametric direction r. 
The construction of T 1 y due to a set of refinement points f of the first level is summarized in 
Algorithm[T] Since the nodes of Ao are represented in homogeneous coordinates with w — 1, the 
transformation to the nodes A| , of local elements T|y can be expressed by a matrix product as 

Ti,j := Aj ;i - = Ty/Ao- (A.3) 

If there are refinement points fg of a higher level, i.e. £ > 1, within an f| / additional transformation 
matrices T 2 ,; are constructed based on A 1 , and fg. The resulting local elements fy., are given by 


hi '■= A 2,i = T 2 ,/Ai j - = T 2 , ; -TijA 0 . 


(A. 4) 


/v v 

The accumulated transformation matrices Ty * relate the final Xg A due to all refinement level to the 
initial knot span Ao- 

T^:=Ay j( - = T> i( -A 0 with T T m „ (A.5) 

m(zL 

where L denotes an index set of all levels defining Tgj which is ordered decreasingly. The set up 
of T g i is described in Algorithm^ 
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Algorithm 1 Set up transformation matrices for the next level 
Require: Node matrix A and related refinement points f of the subsequent level £ 
1: l n = A[l,3] — A[l, 1] 

2: l r2 = A[2,3] — A[2,1] 

3: initiali z e temporary knot vectors A/ and Ay 
4: A/ A[l,k],k =1,3 
5: Ay 4 — A[2, k],k = 1,3 
6: for all f do 
7: A/ G- fi[ 1] 

8: Ay <— T*i [2j 

9: 1 1 = length of each non-zero knot spans of A/ 

10: Zy = length of each non-zero knot spans of Ay 

11: initiali z e array a-/ for transformation matrices T 

12: t,- 2 = 0 

13: for all lj Elj do 

14: zy, = 0 

15: for all /, G h do 

16: T = dmg(li/l n ,lj/lr 2 , 1) 

17: T[1,3]=A[1,1](1 -li/l ri )+t ri 

18: T [2, 3] = A[2, 1] (l — lj/lr 2 ) +tr 2 

19: a T <r- T 

20: t r j = -|- 

21: t/-2 — t r2 + lj 

22: return ay 
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Algorithm 2 Hierarchical refinement 

Require: Node matrix Ao of an element f and refinement points f of all levels 
1: initiali z e array ay for transformation matrices T (. e M 3x3 
2: ay <-T 0 = diag (1,1,1) 

3: for all refinement levels i do 
4: initialize temporary array by for 

5: for all Tg_ l m e ay do [> loop over all computed local elements if | 

1 m ~ T^_ 1)W Ao 

7: initialize temporary array Cf for refinement points ; 

8: for all f(, n e fy do 

9: if fj, : n is inside | m related to A/_ ] m then 

10: Cf <- f Ln 

11: ifcy = 0then 

12: b ? <— T £ _i ;m 

13: else 

14: cy = array of set up by Algorithm [T] with A^_j m and c f 

15: for all T^ r e Cj do 

16: by i T£ r T£_i m 

17: a y — by 

18: return ay 
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